Machine learning project

Authors: Jose Pérez Cano & Álvaro Ribot Barrado

0. Libraries

install.packages("klaR")
install.packages("TunePareto")
install.packages("rgl")
install.packages("glmnet")
install.packages("ca")
# LDA/ QDA
library(MASS)

# RDA
library(klaR)

# Multinomial
library(nnet)

# Cross-Validation
library(TunePareto)

# Naive Bayes
library(e1071)

# k-NN
library(class)

# 3d
library(rgl)

# LASSO
library(Matrix)
library(glmnet)
Loading required package: foreach
Loaded glmnet 2.0-16
# Correspondence analysis
library(ca)

1. Read data

set.seed(2105)
setwd("../data")
The working directory was changed to /Users/joseperezcano/Desktop/CFIS segundo curso/AA1/Project/data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.
clev <- read.csv("cleveland.csv", header=F)
hung <- read.csv("hungarian.csv", header=F)
va <- read.csv("long-beach-va.csv", header=F)
switz <- read.csv("switzerland.csv", header=F)

clev$location <- "cleveland"
hung$location <- "hungarian"
va$location <- "long-beach-va"
switz$location <- "switzerland"

heart1 <- rbind(clev, hung)
heart2 <- rbind(va, switz)
heart <- rbind(heart1, heart2)
head(heart)

2. Preprocess data

We apply clustering and several plotting techniques to have an idea of the dataset. In case there are NAs we will use k-NN for imputation.

To extract new features we will apply PCA and FDA and keep this new features and components apart.

# It says which columns are all missings
# The index are returned in negative to eliminate them
na.columns <- function(dd){
  rmlist <- c()
  for (i in 1:ncol(dd)){
    if (min(dd[,i]) == max(dd[,i]) & min(dd[,i])==-9){
      rmlist <- c(rmlist, i)
    }
  }
  -rmlist
}

clev <- clev[,na.columns(clev)]

# Returns columns with more NA than a given threshold, also in negative
much.na.cols <- function(dd, threshold){
  rmlist <- c()
  for (i in 1:ncol(dd)){
    if (sum(dd[,i]==-9) > threshold){
      rmlist <- c(rmlist, i)
    }
  }   
  -rmlist  
}

clev <- clev[, much.na.cols(clev, 60)]

# Applies k-nearest neighbour imputation for a given variable
knn.imputation = function (dd, variable, varname, k)
{  
  aux = subset (dd, select = names(dd)[names(dd) != varname])
  aux1 = aux[!is.na(variable),]
  aux2 = aux[is.na(variable),]

  # Neither of aux1, aux2 can contain NAs
  knn.inc = knn (aux1,aux2, variable[!is.na(variable)], k)
  variable[is.na(variable)] = knn.inc
  variable
}

# This are the variables which values where substituted by dummy values.
dummy <- c("V1", "V2", "V36", "V69", "V70", "V71", "V72", "V73", "V28", "location")
clev <- clev[,!(names(clev) %in% dummy)]

# knn imputation for clev
na.names <- names(clev)[-much.na.cols(clev, 0)]
for (name in na.names){
  clev[, name][clev[, name] == -9] <- NA
  clev[, name] <- knn.imputation(clev, clev[,name], name, 7)
}

Now, we analyse the correlations among variables since it will have impact on later models.

corr.factors <- cor(clev)
which(abs(corr.factors)-diag(diag(corr.factors))>0.9, arr.ind=T)
    row col
V55  34  12
V57  36  14
V31  21  20
V29  20  21
V20  12  34
V22  14  36
rm.correlated <- c("V57", "V55")
clev <- clev[,!(names(clev) %in% rm.correlated)]

factores <- c("V58", "V4", "V9", "V16", "V18", "V19", "V20", "V21", "V22", "V23", "V24", "V25", "V26", "V27", "V38", "V39", "V41", "V51", "V56", "V11", "V59", "V60", "V61", "V63", "V65", "V67", "V68")
for (f in factores){
  clev[,f] <- as.factor(clev[,f])
}

# Dummy level gets replaced
clev$V25[clev$V25 == 2] <- 1
clev$V25 <- droplevels(clev$V25)
levels(clev$V25)
[1] "0" "1"

This is the dataset after the treatment for missings.

summary(clev)
       V3        V4      V9           V10        V11    
 Min.   :29.00   0: 91   1: 22   Min.   : 94.0   0:108  
 1st Qu.:48.00   1:191   2: 43   1st Qu.:120.0   1:174  
 Median :55.00           3: 84   Median :130.0          
 Mean   :54.41           4:133   Mean   :131.6          
 3rd Qu.:61.00                   3rd Qu.:140.0          
 Max.   :77.00                   Max.   :200.0          
                                                        
      V12             V14             V15        V16     V18    
 Min.   :126.0   Min.   : 0.00   Min.   : 0.00   0:240   0:107  
 1st Qu.:213.0   1st Qu.: 0.00   1st Qu.: 0.00   1: 42   1:175  
 Median :244.0   Median :10.00   Median :15.00                  
 Mean   :249.1   Mean   :16.64   Mean   :15.08                  
 3rd Qu.:277.0   3rd Qu.:30.00   3rd Qu.:30.00                  
 Max.   :564.0   Max.   :99.00   Max.   :54.00                  
                                                                
 V19          V20          V21      V22     V23     V24    
 0:138   1      :38   4      : 13   81:67   0:271   0:186  
 1:  2   12     :33   15     : 13   82:96   1: 11   1: 95  
 2:142   2      :31   20     : 13   83:86           2:  1  
         8      :30   13     : 12   84:33                  
         6      :26   16     : 12                          
         11     :25   21     : 12                          
         (Other):99   (Other):207                          
 V25     V26     V27          V29              V31        
 0:211   0:252   0:248   Min.   : 1.800   Min.   : 3.000  
 1: 71   1: 30   1: 34   1st Qu.: 6.500   1st Qu.: 7.000  
                         Median : 8.500   Median : 9.000  
                         Mean   : 8.418   Mean   : 9.754  
                         3rd Qu.:10.075   3rd Qu.:12.000  
                         Max.   :15.000   Max.   :18.000  
                                                          
      V32             V33              V34       
 Min.   : 71.0   Min.   : 40.00   Min.   : 84.0  
 1st Qu.:133.2   1st Qu.: 65.00   1st Qu.:154.0  
 Median :153.5   Median : 74.00   Median :168.0  
 Mean   :149.8   Mean   : 75.12   Mean   :168.1  
 3rd Qu.:165.8   3rd Qu.: 84.00   3rd Qu.:184.0  
 Max.   :202.0   Max.   :119.00   Max.   :232.0  
                                                 
      V35              V37         V38     V39    
 Min.   : 26.00   Min.   : 50.00   0:190   0:276  
 1st Qu.: 70.00   1st Qu.: 80.00   1: 92   1:  6  
 Median : 80.00   Median : 85.00                  
 Mean   : 78.74   Mean   : 84.95                  
 3rd Qu.: 85.00   3rd Qu.: 90.00                  
 Max.   :120.00   Max.   :110.00                  
                                                  
      V40        V41          V43             V44         V51    
 Min.   :0.000   1:135   Min.   : 24.0   Min.   :0.0000   1:  2  
 1st Qu.:0.000   2:129   1st Qu.: 92.0   1st Qu.:0.0000   3:159  
 Median :0.800   3: 18   Median :118.0   Median :0.0000   6: 14  
 Mean   :1.027           Mean   :123.6   Mean   :0.6702   7:107  
 3rd Qu.:1.600           3rd Qu.:152.8   3rd Qu.:1.0000          
 Max.   :6.200           Max.   :270.0   Max.   :3.0000          
                                                                 
      V56      V58     V59     V60     V61     V63     V65    
 5      : 14   0:157   1:270   1:242   1:224   1:238   1:236  
 21     : 14   1: 50   2: 12   2: 40   2: 58   2: 44   2: 46  
 14     : 12   2: 31                                          
 1      : 11   3: 32                                          
 17     : 11   4: 12                                          
 30     : 11                                                  
 (Other):209                                                  
 V67     V68    
 1:233   1:246  
 2: 49   2: 36  
                
                
                
                
                

2.1 Visualizations

par(mfrow = c(2,3))
for(i in 1:ncol(clev)){
  if (!is.factor(clev[,i])) {
    hist(clev[,i], main = names(clev)[i], xlab="Values")
  }
}

par(mfrow = c(2,3))
for(i in 1:ncol(clev)){
  if (!is.factor(clev[,i])) {
    boxplot(clev[,i], xlab = names(clev)[i])
  }
}

par(mfrow = c(3,3))
for(i in 1:ncol(clev)){
  if (is.factor(clev[,i])) {
    hist(as.numeric(as.character(clev[,i])), main = names(clev)[i], xlab="Values")
  }
}

names_num <- c()
for(i in 1:ncol(clev)){
  if (!is.factor(clev[,i])) {
    names_num <- c(names_num, i)
  }
}
clev_numeric <- clev[,names_num]

clev_cor <- cor(clev_numeric)
which(clev_cor > 0.5 & clev_cor < 1, arr.ind = TRUE)
    row col
V34  10   2
V37  12   2
V15   5   4
V14   4   5
V31   7   6
V29   6   7
V32   8   7
V31   7   8
V10   2  10
V37  12  11
V10   2  12
V35  11  12
which(-clev_cor > 0.5 & -clev_cor < 1, arr.ind = TRUE)
     row col

2.2 Modification of values

par(mfrow = c(2,3))
for(i in 1:length(clev_numeric)){
    qqnorm(clev_numeric[,i], main = c("Q-Q Plot: ", names(clev_numeric)[i]))
    qqline(clev_numeric[,i], col=2)
}

Boxcox

par(mfrow = c(2,3))
for(i in 1:length(clev_numeric)){
  boxcox(lm(clev_numeric[,i]-min(clev_numeric[,i])+1e-6~1),lambda = seq(-1, 1.5, by=0.1), xlab = c(names(clev_numeric))[i])
}

par(mfrow = c(1,3))
#we treat them as special cases
aux <- clev_numeric[,"V14"]
aux <- aux[aux !=0]
boxcox(lm(aux~1),lambda = seq(-2, 1.5, by=0.1))
aux <- clev_numeric[,"V15"]
aux <- aux[aux !=0]
boxcox(lm(aux~1),lambda = seq(-2, 1.5, by=0.1))
aux <- clev_numeric[,"V40"]
aux <- aux[aux !=0]
boxcox(lm(aux~1),lambda = seq(-2, 1.5, by=0.1))

clev_sqrt <- c("V10", "V12", "V31", "V43")
clev_sqrt_especial <- c("V14", "V40")
clev_box <- clev_numeric
#box-cox transformation
for (i in 1:ncol(clev_box)){
  if (names(clev_box)[i] %in% clev_sqrt) {
    clev_box[,i] <- 2*sqrt(clev_box[,i]-min(clev_box[,i])+1e-6)
  } else if (names(clev_box)[i] %in% clev_sqrt_especial){
    clev_box[,i] <- 2*sqrt(clev_box[,i])
  }
}
clev_box <- data.frame(scale(clev_box, scale = T)) #standarized
names(clev_box) <- names(clev_numeric)
clev_box$V14 <- clev_box$V14 - min(clev_box$V14)
clev_box$V15 <- clev_box$V15 - min(clev_box$V15)
clev_box$V40 <- clev_box$V40 - min(clev_box$V40)
par(mfrow = c(2,3))
for(i in 1:ncol(clev_box)){
    qqnorm(clev_box[,i], main = c("Q-Q Plot: ", names(clev_box)[i]))
    qqline(clev_box[,i], col=2)
}

For further models, normalisation will be useful.

And this is the final dataset after processing it.

clev[, names(clev_box)] <- clev_box[,names(clev_box)]
summary(clev)
       V3        V4      V9           V10        V11    
 Min.   :29.00   0: 91   1: 22   Min.   : 94.0   0:108  
 1st Qu.:48.00   1:191   2: 43   1st Qu.:120.0   1:174  
 Median :55.00           3: 84   Median :130.0          
 Mean   :54.41           4:133   Mean   :131.6          
 3rd Qu.:61.00                   3rd Qu.:140.0          
 Max.   :77.00                   Max.   :200.0          
                                                        
      V12             V14             V15        V16     V18    
 Min.   :126.0   Min.   : 0.00   Min.   : 0.00   0:240   0:107  
 1st Qu.:213.0   1st Qu.: 0.00   1st Qu.: 0.00   1: 42   1:175  
 Median :244.0   Median :10.00   Median :15.00                  
 Mean   :249.1   Mean   :16.64   Mean   :15.08                  
 3rd Qu.:277.0   3rd Qu.:30.00   3rd Qu.:30.00                  
 Max.   :564.0   Max.   :99.00   Max.   :54.00                  
                                                                
 V19          V20          V21      V22     V23     V24    
 0:138   1      :38   4      : 13   81:67   0:271   0:186  
 1:  2   12     :33   15     : 13   82:96   1: 11   1: 95  
 2:142   2      :31   20     : 13   83:86           2:  1  
         8      :30   13     : 12   84:33                  
         6      :26   16     : 12                          
         11     :25   21     : 12                          
         (Other):99   (Other):207                          
 V25     V26     V27          V29              V31        
 0:211   0:252   0:248   Min.   : 1.800   Min.   : 3.000  
 1: 71   1: 30   1: 34   1st Qu.: 6.500   1st Qu.: 7.000  
                         Median : 8.500   Median : 9.000  
                         Mean   : 8.418   Mean   : 9.754  
                         3rd Qu.:10.075   3rd Qu.:12.000  
                         Max.   :15.000   Max.   :18.000  
                                                          
      V32             V33              V34       
 Min.   : 71.0   Min.   : 40.00   Min.   : 84.0  
 1st Qu.:133.2   1st Qu.: 65.00   1st Qu.:154.0  
 Median :153.5   Median : 74.00   Median :168.0  
 Mean   :149.8   Mean   : 75.12   Mean   :168.1  
 3rd Qu.:165.8   3rd Qu.: 84.00   3rd Qu.:184.0  
 Max.   :202.0   Max.   :119.00   Max.   :232.0  
                                                 
      V35              V37         V38     V39    
 Min.   : 26.00   Min.   : 50.00   0:190   0:276  
 1st Qu.: 70.00   1st Qu.: 80.00   1: 92   1:  6  
 Median : 80.00   Median : 85.00                  
 Mean   : 78.74   Mean   : 84.95                  
 3rd Qu.: 85.00   3rd Qu.: 90.00                  
 Max.   :120.00   Max.   :110.00                  
                                                  
      V40        V41          V43             V44         V51    
 Min.   :0.000   1:135   Min.   : 24.0   Min.   :0.0000   1:  2  
 1st Qu.:0.000   2:129   1st Qu.: 92.0   1st Qu.:0.0000   3:159  
 Median :0.800   3: 18   Median :118.0   Median :0.0000   6: 14  
 Mean   :1.027           Mean   :123.6   Mean   :0.6702   7:107  
 3rd Qu.:1.600           3rd Qu.:152.8   3rd Qu.:1.0000          
 Max.   :6.200           Max.   :270.0   Max.   :3.0000          
                                                                 
      V56      V58     V59     V60     V61     V63     V65    
 5      : 14   0:157   1:270   1:242   1:224   1:238   1:236  
 21     : 14   1: 50   2: 12   2: 40   2: 58   2: 44   2: 46  
 14     : 12   2: 31                                          
 1      : 11   3: 32                                          
 17     : 11   4: 12                                          
 30     : 11                                                  
 (Other):209                                                  
 V67     V68    
 1:233   1:246  
 2: 49   2: 36  
                
                
                
                
                

2.3. Feature extraction

Separe train and test data, seed for reproducibility.

set.seed(2000)
n <- nrow(clev)
train.lenght <- round(2*n/3)

clev <- clev[sample(n),]
train <- clev[1:train.lenght,]
test <- clev[(train.lenght+1):n,]
names_num <- c()
for(i in 1:ncol(train)){
  if (!is.factor(train[,i])) {
    names_num <- c(names_num, i)
  }
}
train_num <- train[,names_num]

Extract PCA features.

pca <- princomp(train_num)
screeplot(pca)

summary(pca)
Importance of components:
                          Comp.1    Comp.2    Comp.3     Comp.4
Standard deviation     1.8444732 1.5577035 1.3875402 1.21213461
Proportion of Variance 0.2297381 0.1638543 0.1300108 0.09921788
Cumulative Proportion  0.2297381 0.3935924 0.5236032 0.62282104
                           Comp.5     Comp.6     Comp.7    Comp.8
Standard deviation     1.01035405 0.97964763 0.85803949 0.7938423
Proportion of Variance 0.06893431 0.06480791 0.04971676 0.0425556
Cumulative Proportion  0.69175535 0.75656326 0.80628002 0.8488356
                           Comp.9    Comp.10    Comp.11
Standard deviation     0.73133601 0.69242679 0.62465621
Proportion of Variance 0.03611787 0.03237695 0.02634938
Cumulative Proportion  0.88495350 0.91733045 0.94367983
                          Comp.12    Comp.13    Comp.14
Standard deviation     0.57654420 0.50756673 0.43609672
Proportion of Variance 0.02244675 0.01739701 0.01284263
Cumulative Proportion  0.96612658 0.98352359 0.99636622
                           Comp.15
Standard deviation     0.231971913
Proportion of Variance 0.003633784
Cumulative Proportion  1.000000000

Fp <- pca$scores
Gs <- pca$loadings

Fs <- Fp %*% diag(1/pca$sdev)
Gp <- Gs %*% diag(pca$sdev) * 2

col.class <- as.numeric(train$V58)
col.class[col.class==1] <- "red"
col.class[col.class==2] <- "green"
col.class[col.class==3] <- "blue"
col.class[col.class==4] <- "yellow"
col.class[col.class==5] <- "purple"

plot(Fs[,1], Fs[,2], asp=1, col = col.class, xlab = "First principal component", ylab = "Second principal component")
arrows(rep(0,dim(Gs)[1]),rep(0,dim(Gs)[1]), Gp[,1], Gp[,2])
text(Gp[,1], Gp[,2], names(train_num), col = "black")
legend("bottomright", fill=c("red","green", "blue", "yellow", "purple"), legend=c('0','1','2','3', '4'))

V1, V59, V57 creates many problems

problematic <- c("V57", "V59")
train <- train[, !(names(train) %in% problematic)]
fda <- lda(V58~V3+V4+V9+V10+V11+V12+V14+V15+V16+V18+V19+V23+V24+V25+V26+V27+V29+V31+V32+V33+V34+V35+V37+V38+V39+V40+V41+V43+V44+V51+V56+V60+V61+V63+V65+V67+V68, data=train)
#plot(fda)
loadings <- predict(fda)$x
plot(loadings, col = col.class)
legend("bottomright", fill=c("red","green", "blue", "yellow", "purple"), legend=c('0','1','2','3', '4'))

train$LD1 <- loadings[,1]
train$LD2 <- loadings[,2]
train$LD3 <- loadings[,3]
train$LD4 <- loadings[,4]

fda_test <- predict(fda, newdata = test)
test$LD1 <- fda_test$x[,1]
test$LD2 <- fda_test$x[,2]
test$LD3 <- fda_test$x[,3]
test$LD4 <- fda_test$x[,4]

Análisis por correspondencias

ac <- mjca(clev[,names(clev) %in% factores], lambda="Burt")
plot(ac, main="MCA biplot of Burt matrix with data")


ac_ind <- mjca(train[,names(train) %in% factores], lambda="indicator", reti = T)
plot(ac_ind$rowcoord, col = col.class)
legend("bottomright", fill=c("red","green", "blue", "yellow", "purple"), legend=c('0','1','2','3', '4'))


mca.features <- ac_ind$rowcoord

3. Resampling protocol

Principal utility function for doing cross-validation.

First we separate in train and test.

Create the CV function for any model with associated predict and update functions.

cross.validation <- function(data, target, model, times, nfolds, need.class){
  set.seed(0202)
  CV.folds <- generateCVRuns(target, ntimes=times, nfold=nfolds, stratified=TRUE)

  err.total <- c()
  for (i in 1:times){
    err.onetime <- c()
    for (j in 1:nfolds){
      print(paste0("Fold: ", j))
      val <- unlist(CV.folds[[i]][[j]])
      
      tr <- data[-val,]
      va <- data[val,]

      model <- update(model, data=tr)
      pred <- predict(model, newdata=va)
      if (need.class){
        pred <- pred$class
      }
      
      err.table <- table(True=target[val], Predicted=pred)
      err <- 1-sum(diag(err.table))/sum(err.table)
      err.onetime <- c(err.onetime, err)
    }
    err.total <- c(err.total, mean(err.onetime))
    print(paste0("Iteration ", i, ", mean error: ", mean(err.onetime)))
  }
  mean(err.total)
}

Cross-validation for k-Nearest Neighbour

cross.validation.knn <- function(data, target, times, nfolds, K){
  set.seed(0202)
  CV.folds <- generateCVRuns(target, ntimes=times, nfold=nfolds, stratified=TRUE)

  err.total <- c()
  for (i in 1:times){
    err.onetime <- c()
    for (j in 1:nfolds){
      val <- unlist(CV.folds[[i]][[j]])
      
      tr <- data[-val,]
      va <- data[val,]

      pred <- knn(tr, va, target[-val], k = K)

      err.table <- table(True=target[val], Predicted=pred)
      err <- 1-sum(diag(err.table))/sum(err.table)
      err.onetime <- c(err.onetime, err)
    }
    err.total <- c(err.total, mean(err.onetime))
  }
  mean(err.total)
}

Cross-validation Naive-Bayes

cross.validation.naive <- function(data, target, model, times, nfolds){
  set.seed(0202)
  CV.folds <- generateCVRuns(target, ntimes=times, nfold=nfolds, stratified=TRUE)

  err.total <- c()
  for (i in 1:times){
    err.onetime <- c()
    for (j in 1:nfolds){
      print(paste0("Fold: ", j))
      val <- unlist(CV.folds[[i]][[j]])
      
      tr <- data[-val,]
      va <- data[val,]

      model <- naiveBayes(V58~.-LD1-LD2-LD3-LD4, data=tr)
      pred <- predict(model, newdata=va)
      
      err.table <- table(True=target[val], Predicted=pred)
      err <- 1-sum(diag(err.table))/sum(err.table)
      err.onetime <- c(err.onetime, err)
    }
    err.total <- c(err.total, mean(err.onetime))
    print(paste0("Iteration ", i, ", mean error: ", mean(err.onetime)))
  }
  mean(err.total)
}

4. Models

The models we are going to use are: - LDA - QDA - RDA - k-NN - Naïve Bayes - GLM

rda.model <- rda(V58~V3+V4+V9+V10+V11+V12+V14+V15+V16+V18+V19+V20+V21+V22+V23+V24+V25+V26+V27+V29+V31+V32+V33+V34+V35+V37+V38+V39+V40+V41+V43+V44+V51+V56+V60+V61+V63+V65+V67+V68, data=train)
naive.model <- naiveBayes(V58~V3+V4+V9+V10+V11+V12+V14+V15+V16+V18+V19+V20+V21+V22+V23+V24+V25+V26+V27+V29+V31+V32+V33+V34+V35+V37+V38+V39+V40+V41+V43+V44+V51+V56+V60+V61+V63+V65+V67+V68, data=train)
cross.validation(train, train$V58, rda.model, 10, 10, T)
rda.model.fda <- rda(V58~.,data=train)
cross.validation(train, train$V58, rda.model.fda, 10, 10, T)
cross.validation.naive(train, train$V58, naive.model, 10, 10)
err <- c()
for (k in 1:20){
  err <- c(err, cross.validation.knn(train, train$V58, 10,10, k))
}
plot(err, type = "l")

err
 [1] 0.1519883 0.1965205 0.1950585 0.2004678 0.2073977 0.2152924
 [7] 0.2184795 0.2276023 0.2313450 0.2323392 0.2419883 0.2456433
[13] 0.2642398 0.2690351 0.2743275 0.2791228 0.2871053 0.2955556
[19] 0.3003801 0.3050877
cross.validation.knn(train, train$V58, 10, 10, 1)

multinomial.model <- multinom(V58~., data=train)

cross.validation(train, train$V58, multinomial.model, 10, 10, F)

multinomial.model.step <- step(multinomial.model)

cross.validation(train, train$V58, multinomial.model.step, 10, 10, F)

multinomial.model.noFDA <- multinom(V58~.-LD1-LD2-LD3-LD4, data=train)

cross.validation(train, train$V58, multinomial.model.noFDA, 10, 10, F)

multinomial.model.noFDA.step <- step(multinomial.model.noFDA)

cross.validation(train, train$V58, multinomial.model.noFDA.step, 10, 10, F)

Test error

rda.model <- update(rda.model.fda, data=train)
pred.test <- predict(rda.model.fda, test)
pred.test <- pred.test$class
(err.table <- table(True=test$V58, Pred=pred.test))
    Pred
True  0  1  2  3  4
   0 40  5  0  0  0
   1  6 13  1  0  2
   2  0  2  6  4  1
   3  0  0  1  9  0
   4  0  0  0  1  3
(err.test <- 1-sum(diag(err.table))/sum(err.table))
[1] 0.2446809

Parte 2

We will now use model on a dataset not yet investigated: Hungary.

hung <- read.table("../data/processed.hungarian.data", header=F, sep=',', na.strings="?")
summary(hung)
       V1              V2               V3              V4       
 Min.   :28.00   Min.   :0.0000   Min.   :1.000   Min.   : 92.0  
 1st Qu.:42.00   1st Qu.:0.0000   1st Qu.:2.000   1st Qu.:120.0  
 Median :49.00   Median :1.0000   Median :3.000   Median :130.0  
 Mean   :47.83   Mean   :0.7245   Mean   :2.983   Mean   :132.6  
 3rd Qu.:54.00   3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:140.0  
 Max.   :66.00   Max.   :1.0000   Max.   :4.000   Max.   :200.0  
                                                  NA's   :1      
       V5              V6                V7        
 Min.   : 85.0   Min.   :0.00000   Min.   :0.0000  
 1st Qu.:209.0   1st Qu.:0.00000   1st Qu.:0.0000  
 Median :243.0   Median :0.00000   Median :0.0000  
 Mean   :250.8   Mean   :0.06993   Mean   :0.2184  
 3rd Qu.:282.5   3rd Qu.:0.00000   3rd Qu.:0.0000  
 Max.   :603.0   Max.   :1.00000   Max.   :2.0000  
 NA's   :23      NA's   :8         NA's   :1       
       V8              V9              V10        
 Min.   : 82.0   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:122.0   1st Qu.:0.0000   1st Qu.:0.0000  
 Median :140.0   Median :0.0000   Median :0.0000  
 Mean   :139.1   Mean   :0.3038   Mean   :0.5861  
 3rd Qu.:155.0   3rd Qu.:1.0000   3rd Qu.:1.0000  
 Max.   :190.0   Max.   :1.0000   Max.   :5.0000  
 NA's   :1       NA's   :1                        
      V11             V12           V13             V14        
 Min.   :1.000   Min.   :0     Min.   :3.000   Min.   :0.0000  
 1st Qu.:2.000   1st Qu.:0     1st Qu.:5.250   1st Qu.:0.0000  
 Median :2.000   Median :0     Median :6.000   Median :0.0000  
 Mean   :1.894   Mean   :0     Mean   :5.643   Mean   :0.3605  
 3rd Qu.:2.000   3rd Qu.:0     3rd Qu.:7.000   3rd Qu.:1.0000  
 Max.   :3.000   Max.   :0     Max.   :7.000   Max.   :1.0000  
 NA's   :190     NA's   :291   NA's   :266                     

Missing values treatment

rm.var <- c("V11", "V12", "V13")
hung <- hung[,!(names(hung) %in% rm.var)]
hung[is.na(hung)] <- -9
var.imputable <- c("V4", "V5", "V6", "V7", "V8", "V9")
for (nom in var.imputable){
  hung[, nom][hung[, nom] == -9] <- NA
  variable <- hung[, nom]
  hung[, nom] <- knn.imputation(hung, variable, nom, 7)
}
var.factor <- c("V2", "V3", "V6", "V7", "V9", "V14")
for (nom in var.factor){
  hung[, nom] <- as.factor(as.character(hung[,nom]))
}
summary(hung)
       V1        V2      V3            V4              V5       
 Min.   :28.00   0: 81   1: 11   Min.   : 27.0   Min.   :  4.0  
 1st Qu.:42.00   1:213   2:106   1st Qu.:120.0   1st Qu.:198.0  
 Median :49.00           3: 54   Median :130.0   Median :237.0  
 Mean   :47.83           4:123   Mean   :132.2   Mean   :236.6  
 3rd Qu.:54.00                   3rd Qu.:140.0   3rd Qu.:277.0  
 Max.   :66.00                   Max.   :200.0   Max.   :603.0  
 V6      V7            V8        V9           V10         V14    
 0:266   0:235   Min.   : 54.0   0:204   Min.   :0.0000   0:188  
 1: 28   1: 53   1st Qu.:122.0   1: 89   1st Qu.:0.0000   1:106  
         2:  6   Median :140.0   2:  1   Median :0.0000          
                 Mean   :138.8           Mean   :0.5861          
                 3rd Qu.:155.0           3rd Qu.:1.0000          
                 Max.   :190.0           Max.   :5.0000          

Gaussianity treatment

names_num <- c()
for(i in 1:ncol(hung)){
  if (!is.factor(hung[,i])) {
    names_num <- c(names_num, i)
  }
}
hung_numeric <- hung[,names_num]

hung_cor <- cor(hung_numeric)
which(hung_cor > 0.5 & hung_cor < 1, arr.ind = TRUE)
     row col
which(-hung_cor > 0.5 & -hung_cor < 1, arr.ind = TRUE)
     row col
par(mfrow = c(2,3))
for(i in 1:length(hung_numeric)){
    qqnorm(hung_numeric[,i], main = c("Q-Q Plot: ", names(hung_numeric)[i]))
    qqline(hung_numeric[,i], col=2)
}

Boxcox

par(mfrow = c(2,3))
for(i in 1:length(hung_numeric)){
  boxcox(lm(hung_numeric[,i]-min(hung_numeric[,i])+1e-6~1),lambda = seq(-1, 1.5, by=0.1), xlab = c(names(hung_numeric))[i])
}

We see var 10 requires a logaritmic transformation.

qqnorm(hung[,"V10"])
qqline(hung[,"V10"], col=2)

Feature extraction

Separe train and test data, seed for reproducibility.

set.seed(2000)
n <- nrow(hung)
train.lenght <- round(2*n/3)

hung <- hung[sample(n),]
train <- hung[1:train.lenght,]
test <- hung[(train.lenght+1):n,]
names_num <- c()
for(i in 1:ncol(train)){
  if (!is.factor(train[,i])) {
    names_num <- c(names_num, i)
  }
}
train_num <- train[,names_num]

Extract PCA features.

pca <- princomp(train_num)
screeplot(pca)

summary(pca)
Importance of components:
                           Comp.1      Comp.2      Comp.3
Standard deviation     86.0161788 24.47912423 17.80405417
Proportion of Variance  0.8849996  0.07167612  0.03791583
Cumulative Proportion   0.8849996  0.95667568  0.99459151
                            Comp.4       Comp.5
Standard deviation     6.709559731 0.4448637716
Proportion of Variance 0.005384815 0.0000236721
Cumulative Proportion  0.999976328 1.0000000000

Fp <- pca$scores
Gs <- pca$loadings

Fs <- Fp %*% diag(1/pca$sdev)
Gp <- Gs %*% diag(pca$sdev) *0.3

col.class <- as.numeric(train$V14)
col.class[col.class==0] <- "red"
col.class[col.class==1] <- "blue"

plot(Fs[,1], Fs[,2], asp=1, col = col.class, xlab = "First principal component", ylab = "Second principal component")
arrows(rep(0,dim(Gs)[1]),rep(0,dim(Gs)[1]), Gp[,1], Gp[,2])
text(Gp[,1], Gp[,2], names(train_num), col = "black")
legend("bottomright", fill=c("red","blue"), legend=c('0','1'))

fda <- lda(V14~., data=train)
#plot(fda)
loadings <- predict(fda)$x
plot(loadings, col = col.class)
legend("bottomright", fill=c("red","blue"), legend=c('0','1'))

train$LD1 <- loadings[,1]

fda_test <- predict(fda, newdata = test)
test$LD1 <- fda_test$x[,1]

Análisis por correspondencias

ac <- mjca(hung[,names(hung) %in% var.factor], lambda="Burt")
plot(ac, main="MCA biplot of Burt matrix with data")


ac_ind <- mjca(train[,names(train) %in% var.factor], lambda="indicator", reti = T)
plot(ac_ind$rowcoord, col = col.class)
legend("bottomright", fill=c("red","blue"), legend=c('0','1'))


mca.features <- ac_ind$rowcoord

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKCiMgTWFjaGluZSBsZWFybmluZyBwcm9qZWN0CiMjIEF1dGhvcnM6IEpvc2UgUMOpcmV6IENhbm8gJiDDgWx2YXJvIFJpYm90IEJhcnJhZG8KCiMjIyAwLiBMaWJyYXJpZXMKCmBgYHtyfQppbnN0YWxsLnBhY2thZ2VzKCJrbGFSIikKaW5zdGFsbC5wYWNrYWdlcygiVHVuZVBhcmV0byIpCmluc3RhbGwucGFja2FnZXMoInJnbCIpCmluc3RhbGwucGFja2FnZXMoImdsbW5ldCIpCmluc3RhbGwucGFja2FnZXMoImNhIikKYGBgCgpgYGB7cn0KIyBMREEvIFFEQQpsaWJyYXJ5KE1BU1MpCgojIFJEQQpsaWJyYXJ5KGtsYVIpCgojIE11bHRpbm9taWFsCmxpYnJhcnkobm5ldCkKCiMgQ3Jvc3MtVmFsaWRhdGlvbgpsaWJyYXJ5KFR1bmVQYXJldG8pCgojIE5haXZlIEJheWVzCmxpYnJhcnkoZTEwNzEpCgojIGstTk4KbGlicmFyeShjbGFzcykKCiMgM2QKbGlicmFyeShyZ2wpCgojIExBU1NPCmxpYnJhcnkoTWF0cml4KQpsaWJyYXJ5KGdsbW5ldCkKCiMgQ29ycmVzcG9uZGVuY2UgYW5hbHlzaXMKbGlicmFyeShjYSkKYGBgCgoKIyMjIDEuIFJlYWQgZGF0YQoKYGBge3J9CnNldC5zZWVkKDIxMDUpCnNldHdkKCIuLi9kYXRhIikKY2xldiA8LSByZWFkLmNzdigiY2xldmVsYW5kLmNzdiIsIGhlYWRlcj1GKQpodW5nIDwtIHJlYWQuY3N2KCJodW5nYXJpYW4uY3N2IiwgaGVhZGVyPUYpCnZhIDwtIHJlYWQuY3N2KCJsb25nLWJlYWNoLXZhLmNzdiIsIGhlYWRlcj1GKQpzd2l0eiA8LSByZWFkLmNzdigic3dpdHplcmxhbmQuY3N2IiwgaGVhZGVyPUYpCgpjbGV2JGxvY2F0aW9uIDwtICJjbGV2ZWxhbmQiCmh1bmckbG9jYXRpb24gPC0gImh1bmdhcmlhbiIKdmEkbG9jYXRpb24gPC0gImxvbmctYmVhY2gtdmEiCnN3aXR6JGxvY2F0aW9uIDwtICJzd2l0emVybGFuZCIKCmhlYXJ0MSA8LSByYmluZChjbGV2LCBodW5nKQpoZWFydDIgPC0gcmJpbmQodmEsIHN3aXR6KQpoZWFydCA8LSByYmluZChoZWFydDEsIGhlYXJ0MikKaGVhZChoZWFydCkKYGBgCgoKIyMjIDIuIFByZXByb2Nlc3MgZGF0YQpXZSBhcHBseSBjbHVzdGVyaW5nIGFuZCBzZXZlcmFsIHBsb3R0aW5nIHRlY2huaXF1ZXMgdG8gaGF2ZSBhbiBpZGVhIG9mIHRoZSBkYXRhc2V0LiBJbiBjYXNlIHRoZXJlIGFyZSBOQXMgd2Ugd2lsbCB1c2Ugay1OTiBmb3IgaW1wdXRhdGlvbi4gCgpUbyBleHRyYWN0IG5ldyBmZWF0dXJlcyB3ZSB3aWxsIGFwcGx5IFBDQSBhbmQgRkRBIGFuZCBrZWVwIHRoaXMgbmV3IGZlYXR1cmVzIGFuZCBjb21wb25lbnRzIGFwYXJ0LgoKYGBge3J9CiMgSXQgc2F5cyB3aGljaCBjb2x1bW5zIGFyZSBhbGwgbWlzc2luZ3MKIyBUaGUgaW5kZXggYXJlIHJldHVybmVkIGluIG5lZ2F0aXZlIHRvIGVsaW1pbmF0ZSB0aGVtCm5hLmNvbHVtbnMgPC0gZnVuY3Rpb24oZGQpewogIHJtbGlzdCA8LSBjKCkKICBmb3IgKGkgaW4gMTpuY29sKGRkKSl7CiAgICBpZiAobWluKGRkWyxpXSkgPT0gbWF4KGRkWyxpXSkgJiBtaW4oZGRbLGldKT09LTkpewogICAgICBybWxpc3QgPC0gYyhybWxpc3QsIGkpCiAgICB9CiAgfQogIC1ybWxpc3QKfQoKY2xldiA8LSBjbGV2WyxuYS5jb2x1bW5zKGNsZXYpXQoKIyBSZXR1cm5zIGNvbHVtbnMgd2l0aCBtb3JlIE5BIHRoYW4gYSBnaXZlbiB0aHJlc2hvbGQsIGFsc28gaW4gbmVnYXRpdmUKbXVjaC5uYS5jb2xzIDwtIGZ1bmN0aW9uKGRkLCB0aHJlc2hvbGQpewogIHJtbGlzdCA8LSBjKCkKICBmb3IgKGkgaW4gMTpuY29sKGRkKSl7CiAgICBpZiAoc3VtKGRkWyxpXT09LTkpID4gdGhyZXNob2xkKXsKICAgICAgcm1saXN0IDwtIGMocm1saXN0LCBpKQogICAgfQogIH0gICAKICAtcm1saXN0ICAKfQoKY2xldiA8LSBjbGV2WywgbXVjaC5uYS5jb2xzKGNsZXYsIDYwKV0KCiMgQXBwbGllcyBrLW5lYXJlc3QgbmVpZ2hib3VyIGltcHV0YXRpb24gZm9yIGEgZ2l2ZW4gdmFyaWFibGUKa25uLmltcHV0YXRpb24gPSBmdW5jdGlvbiAoZGQsIHZhcmlhYmxlLCB2YXJuYW1lLCBrKQp7ICAKICBhdXggPSBzdWJzZXQgKGRkLCBzZWxlY3QgPSBuYW1lcyhkZClbbmFtZXMoZGQpICE9IHZhcm5hbWVdKQogIGF1eDEgPSBhdXhbIWlzLm5hKHZhcmlhYmxlKSxdCiAgYXV4MiA9IGF1eFtpcy5uYSh2YXJpYWJsZSksXQoKICAjIE5laXRoZXIgb2YgYXV4MSwgYXV4MiBjYW4gY29udGFpbiBOQXMKICBrbm4uaW5jID0ga25uIChhdXgxLGF1eDIsIHZhcmlhYmxlWyFpcy5uYSh2YXJpYWJsZSldLCBrKQogIHZhcmlhYmxlW2lzLm5hKHZhcmlhYmxlKV0gPSBrbm4uaW5jCiAgdmFyaWFibGUKfQoKIyBUaGlzIGFyZSB0aGUgdmFyaWFibGVzIHdoaWNoIHZhbHVlcyB3aGVyZSBzdWJzdGl0dXRlZCBieSBkdW1teSB2YWx1ZXMuCmR1bW15IDwtIGMoIlYxIiwgIlYyIiwgIlYzNiIsICJWNjkiLCAiVjcwIiwgIlY3MSIsICJWNzIiLCAiVjczIiwgIlYyOCIsICJsb2NhdGlvbiIpCmNsZXYgPC0gY2xldlssIShuYW1lcyhjbGV2KSAlaW4lIGR1bW15KV0KCiMga25uIGltcHV0YXRpb24gZm9yIGNsZXYKbmEubmFtZXMgPC0gbmFtZXMoY2xldilbLW11Y2gubmEuY29scyhjbGV2LCAwKV0KZm9yIChuYW1lIGluIG5hLm5hbWVzKXsKICBjbGV2WywgbmFtZV1bY2xldlssIG5hbWVdID09IC05XSA8LSBOQQogIGNsZXZbLCBuYW1lXSA8LSBrbm4uaW1wdXRhdGlvbihjbGV2LCBjbGV2WyxuYW1lXSwgbmFtZSwgNykKfQpgYGAKCk5vdywgd2UgYW5hbHlzZSB0aGUgY29ycmVsYXRpb25zIGFtb25nIHZhcmlhYmxlcyBzaW5jZSBpdCB3aWxsIGhhdmUgaW1wYWN0IG9uIGxhdGVyIG1vZGVscy4KCmBgYHtyfQpjb3JyLmZhY3RvcnMgPC0gY29yKGNsZXYpCndoaWNoKGFicyhjb3JyLmZhY3RvcnMpLWRpYWcoZGlhZyhjb3JyLmZhY3RvcnMpKT4wLjksIGFyci5pbmQ9VCkKCnJtLmNvcnJlbGF0ZWQgPC0gYygiVjU3IiwgIlY1NSIpCmNsZXYgPC0gY2xldlssIShuYW1lcyhjbGV2KSAlaW4lIHJtLmNvcnJlbGF0ZWQpXQoKZmFjdG9yZXMgPC0gYygiVjU4IiwgIlY0IiwgIlY5IiwgIlYxNiIsICJWMTgiLCAiVjE5IiwgIlYyMCIsICJWMjEiLCAiVjIyIiwgIlYyMyIsICJWMjQiLCAiVjI1IiwgIlYyNiIsICJWMjciLCAiVjM4IiwgIlYzOSIsICJWNDEiLCAiVjUxIiwgIlY1NiIsICJWMTEiLCAiVjU5IiwgIlY2MCIsICJWNjEiLCAiVjYzIiwgIlY2NSIsICJWNjciLCAiVjY4IikKZm9yIChmIGluIGZhY3RvcmVzKXsKICBjbGV2WyxmXSA8LSBhcy5mYWN0b3IoY2xldlssZl0pCn0KCiMgRHVtbXkgbGV2ZWwgZ2V0cyByZXBsYWNlZApjbGV2JFYyNVtjbGV2JFYyNSA9PSAyXSA8LSAxCmNsZXYkVjI1IDwtIGRyb3BsZXZlbHMoY2xldiRWMjUpCmxldmVscyhjbGV2JFYyNSkKYGBgCgpUaGlzIGlzIHRoZSBkYXRhc2V0IGFmdGVyIHRoZSB0cmVhdG1lbnQgZm9yIG1pc3NpbmdzLgoKYGBge3J9CnN1bW1hcnkoY2xldikKYGBgCgoKIyMjIyAyLjEgVmlzdWFsaXphdGlvbnMKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bmNvbChjbGV2KSl7CiAgaWYgKCFpcy5mYWN0b3IoY2xldlssaV0pKSB7CiAgICBoaXN0KGNsZXZbLGldLCBtYWluID0gbmFtZXMoY2xldilbaV0sIHhsYWI9IlZhbHVlcyIpCiAgfQp9CmBgYAoKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bmNvbChjbGV2KSl7CiAgaWYgKCFpcy5mYWN0b3IoY2xldlssaV0pKSB7CiAgICBib3hwbG90KGNsZXZbLGldLCB4bGFiID0gbmFtZXMoY2xldilbaV0pCiAgfQp9CmBgYAoKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDMsMykpCmZvcihpIGluIDE6bmNvbChjbGV2KSl7CiAgaWYgKGlzLmZhY3RvcihjbGV2WyxpXSkpIHsKICAgIGhpc3QoYXMubnVtZXJpYyhhcy5jaGFyYWN0ZXIoY2xldlssaV0pKSwgbWFpbiA9IG5hbWVzKGNsZXYpW2ldLCB4bGFiPSJWYWx1ZXMiKQogIH0KfQpgYGAKCgpgYGB7cn0KbmFtZXNfbnVtIDwtIGMoKQpmb3IoaSBpbiAxOm5jb2woY2xldikpewogIGlmICghaXMuZmFjdG9yKGNsZXZbLGldKSkgewogICAgbmFtZXNfbnVtIDwtIGMobmFtZXNfbnVtLCBpKQogIH0KfQpjbGV2X251bWVyaWMgPC0gY2xldlssbmFtZXNfbnVtXQoKY2xldl9jb3IgPC0gY29yKGNsZXZfbnVtZXJpYykKd2hpY2goY2xldl9jb3IgPiAwLjUgJiBjbGV2X2NvciA8IDEsIGFyci5pbmQgPSBUUlVFKQp3aGljaCgtY2xldl9jb3IgPiAwLjUgJiAtY2xldl9jb3IgPCAxLCBhcnIuaW5kID0gVFJVRSkKYGBgCgoKIyMjIyAyLjIgTW9kaWZpY2F0aW9uIG9mIHZhbHVlcwoKYGBge3J9CnBhcihtZnJvdyA9IGMoMiwzKSkKZm9yKGkgaW4gMTpsZW5ndGgoY2xldl9udW1lcmljKSl7CiAgICBxcW5vcm0oY2xldl9udW1lcmljWyxpXSwgbWFpbiA9IGMoIlEtUSBQbG90OiAiLCBuYW1lcyhjbGV2X251bWVyaWMpW2ldKSkKICAgIHFxbGluZShjbGV2X251bWVyaWNbLGldLCBjb2w9MikKfQpgYGAKCgpCb3hjb3gKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bGVuZ3RoKGNsZXZfbnVtZXJpYykpewogIGJveGNveChsbShjbGV2X251bWVyaWNbLGldLW1pbihjbGV2X251bWVyaWNbLGldKSsxZS02fjEpLGxhbWJkYSA9IHNlcSgtMSwgMS41LCBieT0wLjEpLCB4bGFiID0gYyhuYW1lcyhjbGV2X251bWVyaWMpKVtpXSkKfQpgYGAKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDEsMykpCiN3ZSB0cmVhdCB0aGVtIGFzIHNwZWNpYWwgY2FzZXMKYXV4IDwtIGNsZXZfbnVtZXJpY1ssIlYxNCJdCmF1eCA8LSBhdXhbYXV4ICE9MF0KYm94Y294KGxtKGF1eH4xKSxsYW1iZGEgPSBzZXEoLTIsIDEuNSwgYnk9MC4xKSkKYXV4IDwtIGNsZXZfbnVtZXJpY1ssIlYxNSJdCmF1eCA8LSBhdXhbYXV4ICE9MF0KYm94Y294KGxtKGF1eH4xKSxsYW1iZGEgPSBzZXEoLTIsIDEuNSwgYnk9MC4xKSkKYXV4IDwtIGNsZXZfbnVtZXJpY1ssIlY0MCJdCmF1eCA8LSBhdXhbYXV4ICE9MF0KYm94Y294KGxtKGF1eH4xKSxsYW1iZGEgPSBzZXEoLTIsIDEuNSwgYnk9MC4xKSkKYGBgCgoKYGBge3J9CmNsZXZfc3FydCA8LSBjKCJWMTAiLCAiVjEyIiwgIlYzMSIsICJWNDMiKQpjbGV2X3NxcnRfZXNwZWNpYWwgPC0gYygiVjE0IiwgIlY0MCIpCmNsZXZfYm94IDwtIGNsZXZfbnVtZXJpYwojYm94LWNveCB0cmFuc2Zvcm1hdGlvbgpmb3IgKGkgaW4gMTpuY29sKGNsZXZfYm94KSl7CiAgaWYgKG5hbWVzKGNsZXZfYm94KVtpXSAlaW4lIGNsZXZfc3FydCkgewogICAgY2xldl9ib3hbLGldIDwtIDIqc3FydChjbGV2X2JveFssaV0tbWluKGNsZXZfYm94WyxpXSkrMWUtNikKICB9IGVsc2UgaWYgKG5hbWVzKGNsZXZfYm94KVtpXSAlaW4lIGNsZXZfc3FydF9lc3BlY2lhbCl7CiAgICBjbGV2X2JveFssaV0gPC0gMipzcXJ0KGNsZXZfYm94WyxpXSkKICB9Cn0KY2xldl9ib3ggPC0gZGF0YS5mcmFtZShzY2FsZShjbGV2X2JveCwgc2NhbGUgPSBUKSkgI3N0YW5kYXJpemVkCm5hbWVzKGNsZXZfYm94KSA8LSBuYW1lcyhjbGV2X251bWVyaWMpCmNsZXZfYm94JFYxNCA8LSBjbGV2X2JveCRWMTQgLSBtaW4oY2xldl9ib3gkVjE0KQpjbGV2X2JveCRWMTUgPC0gY2xldl9ib3gkVjE1IC0gbWluKGNsZXZfYm94JFYxNSkKY2xldl9ib3gkVjQwIDwtIGNsZXZfYm94JFY0MCAtIG1pbihjbGV2X2JveCRWNDApCmBgYAoKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bmNvbChjbGV2X2JveCkpewogICAgcXFub3JtKGNsZXZfYm94WyxpXSwgbWFpbiA9IGMoIlEtUSBQbG90OiAiLCBuYW1lcyhjbGV2X2JveClbaV0pKQogICAgcXFsaW5lKGNsZXZfYm94WyxpXSwgY29sPTIpCn0KYGBgCgpGb3IgZnVydGhlciBtb2RlbHMsIG5vcm1hbGlzYXRpb24gd2lsbCBiZSB1c2VmdWwuCgpgYGB7cn0KY2xldl9ib3ggPC0gYXBwbHkoY2xldl9ib3gsIDIsIHNjYWxlKQpgYGAKCgpBbmQgdGhpcyBpcyB0aGUgZmluYWwgZGF0YXNldCBhZnRlciBwcm9jZXNzaW5nIGl0LgoKYGBge3J9CmNsZXZbLCBuYW1lcyhjbGV2X2JveCldIDwtIGNsZXZfYm94WyxuYW1lcyhjbGV2X2JveCldCnN1bW1hcnkoY2xldikKYGBgCgoKIyMjIyAyLjMuIEZlYXR1cmUgZXh0cmFjdGlvbgoKU2VwYXJlIHRyYWluIGFuZCB0ZXN0IGRhdGEsIHNlZWQgZm9yIHJlcHJvZHVjaWJpbGl0eS4KCmBgYHtyfQpzZXQuc2VlZCgyMDAwKQpuIDwtIG5yb3coY2xldikKdHJhaW4ubGVuZ2h0IDwtIHJvdW5kKDIqbi8zKQoKY2xldiA8LSBjbGV2W3NhbXBsZShuKSxdCnRyYWluIDwtIGNsZXZbMTp0cmFpbi5sZW5naHQsXQp0ZXN0IDwtIGNsZXZbKHRyYWluLmxlbmdodCsxKTpuLF0KYGBgCgpgYGB7cn0KbmFtZXNfbnVtIDwtIGMoKQpmb3IoaSBpbiAxOm5jb2wodHJhaW4pKXsKICBpZiAoIWlzLmZhY3Rvcih0cmFpblssaV0pKSB7CiAgICBuYW1lc19udW0gPC0gYyhuYW1lc19udW0sIGkpCiAgfQp9CnRyYWluX251bSA8LSB0cmFpblssbmFtZXNfbnVtXQpgYGAKCkV4dHJhY3QgUENBIGZlYXR1cmVzLgoKYGBge3J9CnBjYSA8LSBwcmluY29tcCh0cmFpbl9udW0pCnNjcmVlcGxvdChwY2EpCnN1bW1hcnkocGNhKQpgYGAKCgpgYGB7cn0KYmlwbG90KHBjYSkKYGBgCgpgYGB7cn0KRnAgPC0gcGNhJHNjb3JlcwpHcyA8LSBwY2EkbG9hZGluZ3MKCkZzIDwtIEZwICUqJSBkaWFnKDEvcGNhJHNkZXYpCkdwIDwtIEdzICUqJSBkaWFnKHBjYSRzZGV2KSAqIDIKCmNvbC5jbGFzcyA8LSBhcy5udW1lcmljKHRyYWluJFY1OCkKY29sLmNsYXNzW2NvbC5jbGFzcz09MV0gPC0gInJlZCIKY29sLmNsYXNzW2NvbC5jbGFzcz09Ml0gPC0gImdyZWVuIgpjb2wuY2xhc3NbY29sLmNsYXNzPT0zXSA8LSAiYmx1ZSIKY29sLmNsYXNzW2NvbC5jbGFzcz09NF0gPC0gInllbGxvdyIKY29sLmNsYXNzW2NvbC5jbGFzcz09NV0gPC0gInB1cnBsZSIKCnBsb3QoRnNbLDFdLCBGc1ssMl0sIGFzcD0xLCBjb2wgPSBjb2wuY2xhc3MsIHhsYWIgPSAiRmlyc3QgcHJpbmNpcGFsIGNvbXBvbmVudCIsIHlsYWIgPSAiU2Vjb25kIHByaW5jaXBhbCBjb21wb25lbnQiKQphcnJvd3MocmVwKDAsZGltKEdzKVsxXSkscmVwKDAsZGltKEdzKVsxXSksIEdwWywxXSwgR3BbLDJdKQp0ZXh0KEdwWywxXSwgR3BbLDJdLCBuYW1lcyh0cmFpbl9udW0pLCBjb2wgPSAiYmxhY2siKQpsZWdlbmQoImJvdHRvbXJpZ2h0IiwgZmlsbD1jKCJyZWQiLCJncmVlbiIsICJibHVlIiwgInllbGxvdyIsICJwdXJwbGUiKSwgbGVnZW5kPWMoJzAnLCcxJywnMicsJzMnLCAnNCcpKQpgYGAKCgpWMSwgVjU5LCBWNTcgY3JlYXRlcyBtYW55IHByb2JsZW1zCgpgYGB7cn0KcHJvYmxlbWF0aWMgPC0gYygiVjU3IiwgIlY1OSIpCnRyYWluIDwtIHRyYWluWywgIShuYW1lcyh0cmFpbikgJWluJSBwcm9ibGVtYXRpYyldCmZkYSA8LSBsZGEoVjU4flYzK1Y0K1Y5K1YxMCtWMTErVjEyK1YxNCtWMTUrVjE2K1YxOCtWMTkrVjIzK1YyNCtWMjUrVjI2K1YyNytWMjkrVjMxK1YzMitWMzMrVjM0K1YzNStWMzcrVjM4K1YzOStWNDArVjQxK1Y0MytWNDQrVjUxK1Y1NitWNjArVjYxK1Y2MytWNjUrVjY3K1Y2OCwgZGF0YT10cmFpbikKI3Bsb3QoZmRhKQpsb2FkaW5ncyA8LSBwcmVkaWN0KGZkYSkkeApwbG90KGxvYWRpbmdzLCBjb2wgPSBjb2wuY2xhc3MpCmxlZ2VuZCgiYm90dG9tcmlnaHQiLCBmaWxsPWMoInJlZCIsImdyZWVuIiwgImJsdWUiLCAieWVsbG93IiwgInB1cnBsZSIpLCBsZWdlbmQ9YygnMCcsJzEnLCcyJywnMycsICc0JykpCmBgYAoKCmBgYHtyfQp0cmFpbiRMRDEgPC0gbG9hZGluZ3NbLDFdCnRyYWluJExEMiA8LSBsb2FkaW5nc1ssMl0KdHJhaW4kTEQzIDwtIGxvYWRpbmdzWywzXQp0cmFpbiRMRDQgPC0gbG9hZGluZ3NbLDRdCgpmZGFfdGVzdCA8LSBwcmVkaWN0KGZkYSwgbmV3ZGF0YSA9IHRlc3QpCnRlc3QkTEQxIDwtIGZkYV90ZXN0JHhbLDFdCnRlc3QkTEQyIDwtIGZkYV90ZXN0JHhbLDJdCnRlc3QkTEQzIDwtIGZkYV90ZXN0JHhbLDNdCnRlc3QkTEQ0IDwtIGZkYV90ZXN0JHhbLDRdCmBgYAoKCkFuw6FsaXNpcyBwb3IgY29ycmVzcG9uZGVuY2lhcwoKYGBge3J9CmFjIDwtIG1qY2EoY2xldlssbmFtZXMoY2xldikgJWluJSBmYWN0b3Jlc10sIGxhbWJkYT0iQnVydCIpCnBsb3QoYWMsIG1haW49Ik1DQSBiaXBsb3Qgb2YgQnVydCBtYXRyaXggd2l0aCBkYXRhIikKCmFjX2luZCA8LSBtamNhKHRyYWluWyxuYW1lcyh0cmFpbikgJWluJSBmYWN0b3Jlc10sIGxhbWJkYT0iaW5kaWNhdG9yIiwgcmV0aSA9IFQpCnBsb3QoYWNfaW5kJHJvd2Nvb3JkLCBjb2wgPSBjb2wuY2xhc3MpCmxlZ2VuZCgiYm90dG9tcmlnaHQiLCBmaWxsPWMoInJlZCIsImdyZWVuIiwgImJsdWUiLCAieWVsbG93IiwgInB1cnBsZSIpLCBsZWdlbmQ9YygnMCcsJzEnLCcyJywnMycsICc0JykpCgptY2EuZmVhdHVyZXMgPC0gYWNfaW5kJHJvd2Nvb3JkCmBgYAoKCiMjIyAzLiBSZXNhbXBsaW5nIHByb3RvY29sClByaW5jaXBhbCB1dGlsaXR5IGZ1bmN0aW9uIGZvciBkb2luZyBjcm9zcy12YWxpZGF0aW9uLiAKCkZpcnN0IHdlIHNlcGFyYXRlIGluIHRyYWluIGFuZCB0ZXN0LgoKQ3JlYXRlIHRoZSBDViBmdW5jdGlvbiBmb3IgYW55IG1vZGVsIHdpdGggYXNzb2NpYXRlZCBwcmVkaWN0IGFuZCB1cGRhdGUgZnVuY3Rpb25zLgoKYGBge3J9CmNyb3NzLnZhbGlkYXRpb24gPC0gZnVuY3Rpb24oZGF0YSwgdGFyZ2V0LCBtb2RlbCwgdGltZXMsIG5mb2xkcywgbmVlZC5jbGFzcyl7CiAgc2V0LnNlZWQoMDIwMikKICBDVi5mb2xkcyA8LSBnZW5lcmF0ZUNWUnVucyh0YXJnZXQsIG50aW1lcz10aW1lcywgbmZvbGQ9bmZvbGRzLCBzdHJhdGlmaWVkPVRSVUUpCgogIGVyci50b3RhbCA8LSBjKCkKICBmb3IgKGkgaW4gMTp0aW1lcyl7CiAgICBlcnIub25ldGltZSA8LSBjKCkKICAgIGZvciAoaiBpbiAxOm5mb2xkcyl7CiAgICAgIHByaW50KHBhc3RlMCgiRm9sZDogIiwgaikpCiAgICAgIHZhbCA8LSB1bmxpc3QoQ1YuZm9sZHNbW2ldXVtbal1dKQogICAgICAKICAgICAgdHIgPC0gZGF0YVstdmFsLF0KICAgICAgdmEgPC0gZGF0YVt2YWwsXQoKICAgICAgbW9kZWwgPC0gdXBkYXRlKG1vZGVsLCBkYXRhPXRyKQogICAgICBwcmVkIDwtIHByZWRpY3QobW9kZWwsIG5ld2RhdGE9dmEpCiAgICAgIGlmIChuZWVkLmNsYXNzKXsKICAgICAgICBwcmVkIDwtIHByZWQkY2xhc3MKICAgICAgfQogICAgICAKICAgICAgZXJyLnRhYmxlIDwtIHRhYmxlKFRydWU9dGFyZ2V0W3ZhbF0sIFByZWRpY3RlZD1wcmVkKQogICAgICBlcnIgPC0gMS1zdW0oZGlhZyhlcnIudGFibGUpKS9zdW0oZXJyLnRhYmxlKQogICAgICBlcnIub25ldGltZSA8LSBjKGVyci5vbmV0aW1lLCBlcnIpCiAgICB9CiAgICBlcnIudG90YWwgPC0gYyhlcnIudG90YWwsIG1lYW4oZXJyLm9uZXRpbWUpKQogICAgcHJpbnQocGFzdGUwKCJJdGVyYXRpb24gIiwgaSwgIiwgbWVhbiBlcnJvcjogIiwgbWVhbihlcnIub25ldGltZSkpKQogIH0KICBtZWFuKGVyci50b3RhbCkKfQpgYGAKCgpDcm9zcy12YWxpZGF0aW9uIGZvciBrLU5lYXJlc3QgTmVpZ2hib3VyCgpgYGB7cn0KY3Jvc3MudmFsaWRhdGlvbi5rbm4gPC0gZnVuY3Rpb24oZGF0YSwgdGFyZ2V0LCB0aW1lcywgbmZvbGRzLCBLKXsKICBzZXQuc2VlZCgwMjAyKQogIENWLmZvbGRzIDwtIGdlbmVyYXRlQ1ZSdW5zKHRhcmdldCwgbnRpbWVzPXRpbWVzLCBuZm9sZD1uZm9sZHMsIHN0cmF0aWZpZWQ9VFJVRSkKCiAgZXJyLnRvdGFsIDwtIGMoKQogIGZvciAoaSBpbiAxOnRpbWVzKXsKICAgIGVyci5vbmV0aW1lIDwtIGMoKQogICAgZm9yIChqIGluIDE6bmZvbGRzKXsKICAgICAgdmFsIDwtIHVubGlzdChDVi5mb2xkc1tbaV1dW1tqXV0pCiAgICAgIAogICAgICB0ciA8LSBkYXRhWy12YWwsXQogICAgICB2YSA8LSBkYXRhW3ZhbCxdCgogICAgICBwcmVkIDwtIGtubih0ciwgdmEsIHRhcmdldFstdmFsXSwgayA9IEspCgogICAgICBlcnIudGFibGUgPC0gdGFibGUoVHJ1ZT10YXJnZXRbdmFsXSwgUHJlZGljdGVkPXByZWQpCiAgICAgIGVyciA8LSAxLXN1bShkaWFnKGVyci50YWJsZSkpL3N1bShlcnIudGFibGUpCiAgICAgIGVyci5vbmV0aW1lIDwtIGMoZXJyLm9uZXRpbWUsIGVycikKICAgIH0KICAgIGVyci50b3RhbCA8LSBjKGVyci50b3RhbCwgbWVhbihlcnIub25ldGltZSkpCiAgfQogIG1lYW4oZXJyLnRvdGFsKQp9CmBgYAoKCkNyb3NzLXZhbGlkYXRpb24gTmFpdmUtQmF5ZXMKCmBgYHtyfQpjcm9zcy52YWxpZGF0aW9uLm5haXZlIDwtIGZ1bmN0aW9uKGRhdGEsIHRhcmdldCwgbW9kZWwsIHRpbWVzLCBuZm9sZHMpewogIHNldC5zZWVkKDAyMDIpCiAgQ1YuZm9sZHMgPC0gZ2VuZXJhdGVDVlJ1bnModGFyZ2V0LCBudGltZXM9dGltZXMsIG5mb2xkPW5mb2xkcywgc3RyYXRpZmllZD1UUlVFKQoKICBlcnIudG90YWwgPC0gYygpCiAgZm9yIChpIGluIDE6dGltZXMpewogICAgZXJyLm9uZXRpbWUgPC0gYygpCiAgICBmb3IgKGogaW4gMTpuZm9sZHMpewogICAgICBwcmludChwYXN0ZTAoIkZvbGQ6ICIsIGopKQogICAgICB2YWwgPC0gdW5saXN0KENWLmZvbGRzW1tpXV1bW2pdXSkKICAgICAgCiAgICAgIHRyIDwtIGRhdGFbLXZhbCxdCiAgICAgIHZhIDwtIGRhdGFbdmFsLF0KCiAgICAgIG1vZGVsIDwtIG5haXZlQmF5ZXMoVjU4fi4tTEQxLUxEMi1MRDMtTEQ0LCBkYXRhPXRyKQogICAgICBwcmVkIDwtIHByZWRpY3QobW9kZWwsIG5ld2RhdGE9dmEpCiAgICAgIAogICAgICBlcnIudGFibGUgPC0gdGFibGUoVHJ1ZT10YXJnZXRbdmFsXSwgUHJlZGljdGVkPXByZWQpCiAgICAgIGVyciA8LSAxLXN1bShkaWFnKGVyci50YWJsZSkpL3N1bShlcnIudGFibGUpCiAgICAgIGVyci5vbmV0aW1lIDwtIGMoZXJyLm9uZXRpbWUsIGVycikKICAgIH0KICAgIGVyci50b3RhbCA8LSBjKGVyci50b3RhbCwgbWVhbihlcnIub25ldGltZSkpCiAgICBwcmludChwYXN0ZTAoIkl0ZXJhdGlvbiAiLCBpLCAiLCBtZWFuIGVycm9yOiAiLCBtZWFuKGVyci5vbmV0aW1lKSkpCiAgfQogIG1lYW4oZXJyLnRvdGFsKQp9CmBgYAoKCiMjIyA0LiBNb2RlbHMKVGhlIG1vZGVscyB3ZSBhcmUgZ29pbmcgdG8gdXNlIGFyZTogCiAgLSBMREEKICAtIFFEQQogIC0gUkRBCiAgLSBrLU5OCiAgLSBOYcOvdmUgQmF5ZXMKICAtIEdMTQogCgpgYGB7cn0KcmRhLm1vZGVsIDwtIHJkYShWNTh+VjMrVjQrVjkrVjEwK1YxMStWMTIrVjE0K1YxNStWMTYrVjE4K1YxOStWMjArVjIxK1YyMitWMjMrVjI0K1YyNStWMjYrVjI3K1YyOStWMzErVjMyK1YzMytWMzQrVjM1K1YzNytWMzgrVjM5K1Y0MCtWNDErVjQzK1Y0NCtWNTErVjU2K1Y2MCtWNjErVjYzK1Y2NStWNjcrVjY4LCBkYXRhPXRyYWluKQpuYWl2ZS5tb2RlbCA8LSBuYWl2ZUJheWVzKFY1OH5WMytWNCtWOStWMTArVjExK1YxMitWMTQrVjE1K1YxNitWMTgrVjE5K1YyMCtWMjErVjIyK1YyMytWMjQrVjI1K1YyNitWMjcrVjI5K1YzMStWMzIrVjMzK1YzNCtWMzUrVjM3K1YzOCtWMzkrVjQwK1Y0MStWNDMrVjQ0K1Y1MStWNTYrVjYwK1Y2MStWNjMrVjY1K1Y2NytWNjgsIGRhdGE9dHJhaW4pCmBgYAoKCmBgYHtyfQpjcm9zcy52YWxpZGF0aW9uKHRyYWluLCB0cmFpbiRWNTgsIHJkYS5tb2RlbCwgMTAsIDEwLCBUKQpgYGAKCgpgYGB7cn0KcmRhLm1vZGVsLmZkYSA8LSByZGEoVjU4fi4sZGF0YT10cmFpbikKYGBgCgoKYGBge3J9CmNyb3NzLnZhbGlkYXRpb24odHJhaW4sIHRyYWluJFY1OCwgcmRhLm1vZGVsLmZkYSwgMTAsIDEwLCBUKQpgYGAKCgpgYGB7cn0KY3Jvc3MudmFsaWRhdGlvbi5uYWl2ZSh0cmFpbiwgdHJhaW4kVjU4LCBuYWl2ZS5tb2RlbCwgMTAsIDEwKQpgYGAKCgpgYGB7cn0KZXJyIDwtIGMoKQpmb3IgKGsgaW4gMToyMCl7CiAgZXJyIDwtIGMoZXJyLCBjcm9zcy52YWxpZGF0aW9uLmtubih0cmFpbiwgdHJhaW4kVjU4LCAxMCwxMCwgaykpCn0KYGBgCgoKYGBge3J9CnBsb3QoZXJyLCB0eXBlID0gImwiKQplcnIKYGBgCgoKYGBge3J9CmNyb3NzLnZhbGlkYXRpb24ua25uKHRyYWluLCB0cmFpbiRWNTgsIDEwLCAxMCwgMSkKCm11bHRpbm9taWFsLm1vZGVsIDwtIG11bHRpbm9tKFY1OH4uLCBkYXRhPXRyYWluKQoKY3Jvc3MudmFsaWRhdGlvbih0cmFpbiwgdHJhaW4kVjU4LCBtdWx0aW5vbWlhbC5tb2RlbCwgMTAsIDEwLCBGKQoKbXVsdGlub21pYWwubW9kZWwuc3RlcCA8LSBzdGVwKG11bHRpbm9taWFsLm1vZGVsKQoKY3Jvc3MudmFsaWRhdGlvbih0cmFpbiwgdHJhaW4kVjU4LCBtdWx0aW5vbWlhbC5tb2RlbC5zdGVwLCAxMCwgMTAsIEYpCgptdWx0aW5vbWlhbC5tb2RlbC5ub0ZEQSA8LSBtdWx0aW5vbShWNTh+Li1MRDEtTEQyLUxEMy1MRDQsIGRhdGE9dHJhaW4pCgpjcm9zcy52YWxpZGF0aW9uKHRyYWluLCB0cmFpbiRWNTgsIG11bHRpbm9taWFsLm1vZGVsLm5vRkRBLCAxMCwgMTAsIEYpCgptdWx0aW5vbWlhbC5tb2RlbC5ub0ZEQS5zdGVwIDwtIHN0ZXAobXVsdGlub21pYWwubW9kZWwubm9GREEpCgpjcm9zcy52YWxpZGF0aW9uKHRyYWluLCB0cmFpbiRWNTgsIG11bHRpbm9taWFsLm1vZGVsLm5vRkRBLnN0ZXAsIDEwLCAxMCwgRikKYGBgCgoKIyMgVGVzdCBlcnJvcgoKYGBge3J9CnJkYS5tb2RlbCA8LSB1cGRhdGUocmRhLm1vZGVsLmZkYSwgZGF0YT10cmFpbikKcHJlZC50ZXN0IDwtIHByZWRpY3QocmRhLm1vZGVsLmZkYSwgdGVzdCkKcHJlZC50ZXN0IDwtIHByZWQudGVzdCRjbGFzcwooZXJyLnRhYmxlIDwtIHRhYmxlKFRydWU9dGVzdCRWNTgsIFByZWQ9cHJlZC50ZXN0KSkKKGVyci50ZXN0IDwtIDEtc3VtKGRpYWcoZXJyLnRhYmxlKSkvc3VtKGVyci50YWJsZSkpCmBgYAoKIyBQYXJ0ZSAyCgpXZSB3aWxsIG5vdyB1c2UgbW9kZWwgb24gYSBkYXRhc2V0IG5vdCB5ZXQgaW52ZXN0aWdhdGVkOiBIdW5nYXJ5LgoKYGBge3J9Cmh1bmcgPC0gcmVhZC50YWJsZSgiLi4vZGF0YS9wcm9jZXNzZWQuaHVuZ2FyaWFuLmRhdGEiLCBoZWFkZXI9Riwgc2VwPScsJywgbmEuc3RyaW5ncz0iPyIpCnN1bW1hcnkoaHVuZykKYGBgCgojIyBNaXNzaW5nIHZhbHVlcyB0cmVhdG1lbnQKCmBgYHtyfQpybS52YXIgPC0gYygiVjExIiwgIlYxMiIsICJWMTMiKQpodW5nIDwtIGh1bmdbLCEobmFtZXMoaHVuZykgJWluJSBybS52YXIpXQpgYGAKCmBgYHtyfQpodW5nW2lzLm5hKGh1bmcpXSA8LSAtOQp2YXIuaW1wdXRhYmxlIDwtIGMoIlY0IiwgIlY1IiwgIlY2IiwgIlY3IiwgIlY4IiwgIlY5IikKZm9yIChub20gaW4gdmFyLmltcHV0YWJsZSl7CiAgaHVuZ1ssIG5vbV1baHVuZ1ssIG5vbV0gPT0gLTldIDwtIE5BCiAgdmFyaWFibGUgPC0gaHVuZ1ssIG5vbV0KICBodW5nWywgbm9tXSA8LSBrbm4uaW1wdXRhdGlvbihodW5nLCB2YXJpYWJsZSwgbm9tLCA3KQp9CmBgYAoKYGBge3J9CnZhci5mYWN0b3IgPC0gYygiVjIiLCAiVjMiLCAiVjYiLCAiVjciLCAiVjkiLCAiVjE0IikKZm9yIChub20gaW4gdmFyLmZhY3Rvcil7CiAgaHVuZ1ssIG5vbV0gPC0gYXMuZmFjdG9yKGFzLmNoYXJhY3RlcihodW5nWyxub21dKSkKfQoKIyBMZXZlbCAyIGhhcyBmZXcgdmFsdWVzCmh1bmckVjdbaHVuZyRWNyA9PSAyXSA8LSAxCmh1bmckVjcgPC0gZHJvcGxldmVscyhodW5nJFY3KQpsZXZlbHMoaHVuZyRWNykKCiMgTGV2ZWwgMiBoYXMgb25lIGFub21hbCB2YWx1ZQpodW5nJFY5W2h1bmckVjkgPT0gMl0gPC0gMQpodW5nJFY5IDwtIGRyb3BsZXZlbHMoaHVuZyRWOSkKbGV2ZWxzKGh1bmckVjkpCmBgYAoKCmBgYHtyfQpzdW1tYXJ5KGh1bmcpCmBgYAoKIyMgR2F1c3NpYW5pdHkgdHJlYXRtZW50CgpgYGB7cn0KbmFtZXNfbnVtIDwtIGMoKQpmb3IoaSBpbiAxOm5jb2woaHVuZykpewogIGlmICghaXMuZmFjdG9yKGh1bmdbLGldKSkgewogICAgbmFtZXNfbnVtIDwtIGMobmFtZXNfbnVtLCBpKQogIH0KfQpodW5nX251bWVyaWMgPC0gaHVuZ1ssbmFtZXNfbnVtXQoKaHVuZ19jb3IgPC0gY29yKGh1bmdfbnVtZXJpYykKd2hpY2goaHVuZ19jb3IgPiAwLjUgJiBodW5nX2NvciA8IDEsIGFyci5pbmQgPSBUUlVFKQp3aGljaCgtaHVuZ19jb3IgPiAwLjUgJiAtaHVuZ19jb3IgPCAxLCBhcnIuaW5kID0gVFJVRSkKYGBgCgoKYGBge3J9CnBhcihtZnJvdyA9IGMoMiwzKSkKZm9yKGkgaW4gMTpsZW5ndGgoaHVuZ19udW1lcmljKSl7CiAgICBxcW5vcm0oaHVuZ19udW1lcmljWyxpXSwgbWFpbiA9IGMoIlEtUSBQbG90OiAiLCBuYW1lcyhodW5nX251bWVyaWMpW2ldKSkKICAgIHFxbGluZShodW5nX251bWVyaWNbLGldLCBjb2w9MikKfQpgYGAKCgpCb3hjb3gKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bGVuZ3RoKGh1bmdfbnVtZXJpYykpewogIGJveGNveChsbShodW5nX251bWVyaWNbLGldLW1pbihodW5nX251bWVyaWNbLGldKSsxZS02fjEpLGxhbWJkYSA9IHNlcSgtMSwgMS41LCBieT0wLjEpLCB4bGFiID0gYyhuYW1lcyhodW5nX251bWVyaWMpKVtpXSkKfQpgYGAKCldlIHNlZSB2YXIgMTAgcmVxdWlyZXMgYSBsb2dhcml0bWljIHRyYW5zZm9ybWF0aW9uLgoKYGBge3J9Cmh1bmdbLCJWMTAiXSA8LSBsb2coMStodW5nWywiVjEwIl0pCmBgYCAKCmBgYHtyfQpxcW5vcm0oaHVuZ1ssIlYxMCJdKQpxcWxpbmUoaHVuZ1ssIlYxMCJdLCBjb2w9MikKYGBgCgojIyBGZWF0dXJlIGV4dHJhY3Rpb24KClNlcGFyZSB0cmFpbiBhbmQgdGVzdCBkYXRhLCBzZWVkIGZvciByZXByb2R1Y2liaWxpdHkuCgpgYGB7cn0Kc2V0LnNlZWQoMjAwMCkKbiA8LSBucm93KGh1bmcpCnRyYWluLmxlbmdodCA8LSByb3VuZCgyKm4vMykKCmh1bmcgPC0gaHVuZ1tzYW1wbGUobiksXQp0cmFpbiA8LSBodW5nWzE6dHJhaW4ubGVuZ2h0LF0KdGVzdCA8LSBodW5nWyh0cmFpbi5sZW5naHQrMSk6bixdCmBgYAoKYGBge3J9Cm5hbWVzX251bSA8LSBjKCkKZm9yKGkgaW4gMTpuY29sKHRyYWluKSl7CiAgaWYgKCFpcy5mYWN0b3IodHJhaW5bLGldKSkgewogICAgbmFtZXNfbnVtIDwtIGMobmFtZXNfbnVtLCBpKQogIH0KfQp0cmFpbl9udW0gPC0gdHJhaW5bLG5hbWVzX251bV0KYGBgCgpFeHRyYWN0IFBDQSBmZWF0dXJlcy4KCmBgYHtyfQpwY2EgPC0gcHJpbmNvbXAodHJhaW5fbnVtKQpzY3JlZXBsb3QocGNhKQpzdW1tYXJ5KHBjYSkKYGBgCgoKYGBge3J9CmJpcGxvdChwY2EpCmBgYAoKYGBge3J9CkZwIDwtIHBjYSRzY29yZXMKR3MgPC0gcGNhJGxvYWRpbmdzCgpGcyA8LSBGcCAlKiUgZGlhZygxL3BjYSRzZGV2KQpHcCA8LSBHcyAlKiUgZGlhZyhwY2Ekc2RldikgKjAuMwoKY29sLmNsYXNzIDwtIGFzLm51bWVyaWModHJhaW4kVjE0KQpjb2wuY2xhc3NbY29sLmNsYXNzPT0wXSA8LSAicmVkIgpjb2wuY2xhc3NbY29sLmNsYXNzPT0xXSA8LSAiYmx1ZSIKCnBsb3QoRnNbLDFdLCBGc1ssMl0sIGFzcD0xLCBjb2wgPSBjb2wuY2xhc3MsIHhsYWIgPSAiRmlyc3QgcHJpbmNpcGFsIGNvbXBvbmVudCIsIHlsYWIgPSAiU2Vjb25kIHByaW5jaXBhbCBjb21wb25lbnQiKQphcnJvd3MocmVwKDAsZGltKEdzKVsxXSkscmVwKDAsZGltKEdzKVsxXSksIEdwWywxXSwgR3BbLDJdKQp0ZXh0KEdwWywxXSwgR3BbLDJdLCBuYW1lcyh0cmFpbl9udW0pLCBjb2wgPSAiYmxhY2siKQpsZWdlbmQoImJvdHRvbXJpZ2h0IiwgZmlsbD1jKCJyZWQiLCJibHVlIiksIGxlZ2VuZD1jKCcwJywnMScpKQpgYGAKCmBgYHtyfQpmZGEgPC0gbGRhKFYxNH4uLCBkYXRhPXRyYWluKQojcGxvdChmZGEpCmxvYWRpbmdzIDwtIHByZWRpY3QoZmRhKSR4CnBsb3QobG9hZGluZ3MsIGNvbCA9IGNvbC5jbGFzcykKbGVnZW5kKCJib3R0b21yaWdodCIsIGZpbGw9YygicmVkIiwiYmx1ZSIpLCBsZWdlbmQ9YygnMCcsJzEnKSkKYGBgCgoKYGBge3J9CnRyYWluJExEMSA8LSBsb2FkaW5nc1ssMV0KCmZkYV90ZXN0IDwtIHByZWRpY3QoZmRhLCBuZXdkYXRhID0gdGVzdCkKdGVzdCRMRDEgPC0gZmRhX3Rlc3QkeFssMV0KYGBgCgoKQW7DoWxpc2lzIHBvciBjb3JyZXNwb25kZW5jaWFzCgpgYGB7cn0KYWMgPC0gbWpjYShodW5nWyxuYW1lcyhodW5nKSAlaW4lIHZhci5mYWN0b3JdLCBsYW1iZGE9IkJ1cnQiKQpwbG90KGFjLCBtYWluPSJNQ0EgYmlwbG90IG9mIEJ1cnQgbWF0cml4IHdpdGggZGF0YSIpCgphY19pbmQgPC0gbWpjYSh0cmFpblssbmFtZXModHJhaW4pICVpbiUgdmFyLmZhY3Rvcl0sIGxhbWJkYT0iaW5kaWNhdG9yIiwgcmV0aSA9IFQpCnBsb3QoYWNfaW5kJHJvd2Nvb3JkLCBjb2wgPSBjb2wuY2xhc3MpCmxlZ2VuZCgiYm90dG9tcmlnaHQiLCBmaWxsPWMoInJlZCIsImJsdWUiKSwgbGVnZW5kPWMoJzAnLCcxJykpCgptY2EuZmVhdHVyZXMgPC0gYWNfaW5kJHJvd2Nvb3JkCmBgYAoKIyMgCg==